Record ages of non-Markovian scale-invariant random walks

How long is needed for an observable to exceed its previous highest value and establish a new record? This time, known as the age of a record plays a crucial role in quantifying record statistics. Until now, general methods for determining record age statistics have been limited to observations of either independent random variables or successive positions of a Markovian (memoryless) random walk. Here we develop a theoretical framework to determine record age statistics in the presence of memory effects for continuous non-smooth processes that are asymptotically scale-invariant. Our theoretical predictions are confirmed by numerical simulations and experimental realisations of diverse representative non-Markovian random walk models and real time series with memory effects, in fields as diverse as genomics, climatology, hydrology, geology and computer science. Our results reveal the crucial role of the number of records already achieved in time series and change our view on analysing record statistics.

We assume that the random walk (RW) (X t ) t=0,... converges at large time to a continuous and non-smooth scaleinvariant process.Under these conditions, the process is characterized by a walk dimension d w > 1, such that X t ∝ t 1/dw , and the random variable X t /t 1/dw is independent of t.We more generally assume that X t has scaleinvariant increments with a potential ageing, meaning that, for 1 ≪ t ≪ T , X t+T − X T ∝ t 1/d 0 w T α/2 .This defines the aging exponent α [S1, S2] and an effective walk dimension at short times d 0 w ≡ (d −1 w − α/2) −1 .To describe the properties of the RWs, we will use their continuous limit, which requires to introduce a microscopic cut-off either in time or space (as done in [S3]).Relying here on a spatial cut-off ∆x, we consider the time T x0 to first reach the level x 0 = n∆x starting from the origin, as the continuous counterpart of the time T n to break n records.
B. Scale-invariance of the time Tn to break n records and of its increments First, we check that the temporal scale-invariance of X t leads to the spatial scale-invariance of T x0 .Indeed, since X t /t 1/dw is asymptotically independent of time t, X Tx 0 /T 1/dw x0 = x 0 /T 1/dw x0 is asymptotically independent of x 0 (X Tx 0 = x 0 because of continuity).Consequently, T x0 is scale-invariant in the sense that T x0 /x dw 0 is independent of x 0 .Taking ∆x = 1 and x 0 = n∆x = n results in the following form for the cumulative distribution: ) where the scaling function Φ is independent of T and n (see Supplementary Figure 2 for numerical check of the number of records n scale-invariance with time T for representative non-Markovian RWs).
Second, we show that the scale-invariance of the increments X t+T − X T , X t+T − X T ∝ t 1/dw−α/2 T α/2 = t 1/d 0 w T α/2 , implies the scale-invariance of T x0+x1 − T x0 .Consider the random variable X T +t −X T T α/2 t 1/d 0 w , which is independent of t and T as long as 1 ≪ t ≪ T .By replacing T by T x0 and t by T x0+x1 − T x0 , we find that , is asymptotically independent of x 0 and x 1 for 1 ≪ x 1 ≪ x 0 .Using that T x0 /x dw 0 is independent of x 0 , we finally obtain that is independent of x 0 and x 1 .In other words, we have the scale invariance of . It means that the cumulative distribution in the limit 1 ≪ n ≪ m and 1 ≪ T ≪ m dw can be written as where the scaling function Ψ is independent of T , n, and m.This scale-invariance of the time increments is systematically checked numerically for a number of representative non-Markovian RWs in Supplementary Figure 3.

C. Characteristic exponents of the record age distribution
We derive the exponents governing the algebraic decay of the record age distribution S(n, τ ) by elaborating the arguments sketched in the main text.
We note by τ k the k th record age or, in the continuum setting, the time to reach level (k + 1)∆x starting from the first arrival at level k∆x.T n being the time to first reach level n∆x, it is given by the sum of the {τ k }, We make the self-consistent assumption that the tail distribution of τ k is algebraic of exponent smaller than 1.As a consequence, the sum ( S3) is controlled by the largest We now assume (see the next subsection for an analytical justification, as well as Supplementary Figure 4 and Supplementary Figure 5 for systematic numerical checks for a number of representative non-Markovian RWs) that the record ages τ k involved in Eq. ( S4) are asymptotically (n ≫ 1) effectively independent, which leads to We search now the exponents y i and ϵ i characterizing the tail distribution of τ k , S(k, T ) ∝ k −1+ϵ1 T −y1 for T ≪ k dw (regime 1) and S(k, T ) ∝ k −1+ϵ2 T −y2 for T ≫ k dw (regime 2), see Eq. (1) of the main text.
We start with T ≫ n dw , so that all S(k, T ) are in the same regime T ≫ k dw .Using Eq. (S5), The scale-invariance of T n (Eq.( S1)) implies that ϵ 2 /y 2 = d w .
Next, for T ≪ n dw , S(n, T ) ∝ n −1+ϵ1 T −y1 .By splitting the sum over k between the regions k dw ≪ T and k dw ≫ T , we have ) This result, together with Eq. ( S1), yields ϵ 1 /y 1 = d w = ϵ 2 /y 2 .In particular, the survival probability of the record-ages admits a scaling form in T and k: To proceed further, we look at the cumulative distribution of the increments T m+n − T m in the limit 1 ≪ n ≪ m and 1 ≪ T ≪ m dw , Using finally Eq. (S2), we obtain y 1 = 1/d 0 w and 1−ϵ1 y1 = d 0 w − d w , leading to ϵ 1 = d w /d 0 w .This provides the exponents y 1 and ϵ 1 of regime 1 (τ n ≪ n dw ).We note that the algebraic decay of S(n, τ ) holds after the time τ n such that ) This gives the lower bound τ n ≫ n dw−d 0 w in Eq. (1) of the main text for regime 1.
For times τ n ≫ n dw (regime 2), the record age τ n is much larger than the typical time needed to break all the previous records.At this time scale, the memory of the n broken records no longer affects the algebraic time decay of S(n, τ ), which is thus given by the usual persistence exponent θ = y 2 (defined as P(T ≥ τ ) ∝ τ −θ where T is the time needed to reach a given value for the first time).Knowing that ϵ 2 /y 2 = d w , we obtain ϵ 2 = d w θ.
By combining the results derived in this section, we finally obtain Eq. (1) of the main text.
D. General scaling criteria showing that the correlations between record ages are asymptotically irrelevant for their maximum

Effective independence criteria
Here we show that in the calculation of the distribution of the maximum of the record ages, M n ≡ max(τ 0 , . . ., τ n−1 ), the random variables τ k can be treated as effectively independent.To do so, we extend the criteria of Ref. [S5].Dividing the set (τ 0 , . . ., τ n−1 ) in two subsets, (τ 0 , . . ., τ n/2−1 ) and (τ n/2 , . . ., τ n ), this criteria states that the correlations between the τ k are irrelevant if (a) the typical mean cross-correlation between the subsystems is much smaller than (b) the variance of the maximum of record ages of the whole system.
Supplementary Figure 1.Non-smoothness versus smoothness.We compare (a) a non-smooth process to (b) a smooth process by zooming on the trajectory that reaches level n for the first time.For the non-smooth process the time needed to cross level n + ∆x1 (with a ∆x1 ≪ n small enough) starting from level n is significantly smaller than for the smooth process.
Here, however, the random variables {τ k } have diverging second moments, see Eq. ( 1) of the main text.To overcome this difficulty, we extend the criteria of Ref. [S5] which uses variances of random variables and consider fractional moments of {τ k } that converge, see below.Also, the original work of Ref. [S5] considers random variables that are Gaussian and identically distributed.Here we extend the approach to the distribution given in Eq. ( 1) of the main text.

General form of P(τ k ≥ T )
We use dimensional analysis to obtain the scaling form of the tail distribution of record ages, without assuming their independence.There are no additional hypotheses on the RW processes beyond those already requested: (i) continuity, (ii) non-smoothness, and (iii) asymptotic scale-invariance.
We recall that the probability to reach level (k + 1)∆x = n + ∆x starting from k∆x = n at a time larger than T is given by a function of microscopic cut-off ∆x (length), the level number n (length) and T (time).Based on dimensional analysis and scale-invariance of the process, the tail distribution of the record ages takes the functional form P(τ k ≥ T ) = F (∆x/n, T /n dw ).We are interested in the limit ∆x/n ≪ 1 and T ∼ n dw , and consider (without restriction of generality) the limit behaviour with β an exponent that we determine in the following.
For continuous non-smooth processes, when reaching the nth level, the RW crosses the level n infinitely often [S6], see Supplementary Figure 1a for an illustration.In particular, this means that when ∆x is small compared to n, there is a time, much smaller than T ∼ n dw , at which the trajectory crosses level n + ∆x.Thus, the probability to reach level n + ∆x at a time larger than T goes to 0 when ∆x goes to 0, which shows that the exponent β in Eq. ( S11) is strictly positive, β > 0.

Cross-correlations
We wish to find an upper bound of the typical correlations between different τ k of the set (τ 0 , . . ., τ n ), for example, typically quantified by the correlation between τ n/4 and τ 3n/4 .However, computing directly their covariance Cov τ n/4 , τ 3n/4 would lead to a diverging result, because of a heavy-tailed τ k distribution.Therefore, we will make use of the fractional powers of record ages, {τ q k }, where the value of q allows converging second moments of these new variables.

Variance of max
The fluctuations of the maximum of the record ages M n ≡ max(τ 0 , . . ., τ n−1 ) are estimated by relying on the variance of M q n ≡ max(τ q 0 , . . ., τ q n−1 ).Based on Eqs. ( S1) and ( S3) and on the self-consistently checked Eq. ( S4), the maximum's moments (k = 1, 2) can be written as: We choose q such that the integral in Eq. ( S15) converges.Because M n is scale invariant (M n /n dw is a nondeterministic n independent random variable), we finally have

Conclusion
We now can compare the typical cross-correlations to the fluctuations of the maximum, Cov τ q n/4 , τ q 3n/4 ≤ Var(τ q n/4 )Var(τ We conclude that the fluctuations of the maximum of record ages dominate the correlations between two record ages, so that the random variables τ k can be considered as effectively independent. We note that, a priori, the correlations between {τ k } are not negligible if one of the three hypotheses ((i) continuity, (ii) non-smoothness, and (iii) scale-invariance) breaks.While the absence of scale-invariance and of continuity would immediately invalidate calculations of this section, it is not evident for smooth processes.The necessity of the assumption of non smoothness originates from the following argument.When a realisation of a smooth process just reaches the level n, in the following the trajectory can go back to position x < n without crossing level n, see Supplementary Figure 1b.Thus, for smooth processes, the probability for a RW trajectory to reach level n + ∆x at a time T for ∆x arbitrary small is finite, implying β = 0 in Eq. ( S11), which invalidates the estimations after it.

S2. NON-MARKOVIAN RANDOM WALKS (RWS)
A. Definition of the non-Markovian RW models In this subsection, we present the non-Markovian random walk (RW) processes, which are used in Fig. 2 of the main text.These processes encompass the three classes of different statistical statistical mechanisms giving rise to a non-Markovian evolution discussed in the main text and depicted in Fig. 1 of the main text.Supplementary Table 1 provides a summary of their characteristic parameters, namely d w , α, d 0 w and θ.
Supplementary Table 1.Summary of the non-Markovian models considered in this study and of their characteristic parameters.
(a) Fractional Brownian motion (fBm).The fBm is a non-Markovian Gaussian process, with stationary increments.Thus, an fBm X t of Hurst index H is defined by its covariance The steps η t = X t − X t−1 are called fractional Gaussian noise (fGn).Nowadays, the fBm is broadly spread and its implementations could be found in standard packages of python or Wolfram Mathematica.Besides, the survival probability of fBm is characterized by the persistence exponent θ = 1 − H, which was derived in [S8-S10].
(b) Quenched fBm (qfBm).This process is an extension of fBm to quenched initial conditions, which results in non-stationary increment statistics.In particular, it describes the height fluctuations under Gaussian noise of an initially flat interface.Then X t corresponds to the height of the interface at position x = 0, X t = h(0, t), h(x, t) following the Stochastic Differential Equation (SDE) Here η(x, t) is a Gaussian noise with possible spatial correlations.We solve numerically this SDE with a spatial discretization ∆x = 1 and a time discretization ∆t = 0.1.The system is initially flat, h(x, t = 0) = h 0 .The model at z = 2 with space-independent noise is a non-stationary fBm of Hurst exponent H = (1 − 1/z)/2 = 1/4 which corresponds to the continuous limit of a solid-on-solid model [S11].The persistence exponent θ was numerically estimated to be θ = 1.55 ± 0.02 [S7].
(c) Elephant RW (eRW).This process is representative of interactions with its own trajectory.At time t, the step η t is drawn uniformly among all the previous steps η i (i < t) and is reversed with probability β.The persistence exponent was determined in [S12] to be θ = 3/2 − 2β.
(d) Self-attractive walk (SATW).This model is a prototypical example of self-interacting RWs.In the SATW model [S12-S15], the RW at position i jumps to a neighbouring site j = i±1 with probability depending on the number of times n j it has visited site j, where H(0) = 0, H(n > 0) = 1 and β > 0. It was shown in [S12] that the persistence exponent of this RW is given by e −β /2.
(e-f) Exponential self-repelling RW.This is another example of self-interacting RW.In this model, the RW at position i jumps to a neighbouring site j = i ± 1 with probability depending on the number of times n j it has visited site j, where κ and β are two positive real numbers.It was shown [S16] that the walk dimension of such walks is given by d w = 1+κ 2+κ .In the case κ = 1 (the True Self-Avoiding Walk, TSAW, [S17-S19]), the persistence exponent is given by [S15] θ = 1/3, while for κ < 1 (the Sub-Exponential Self-repelling Walk, SESRW, [S16, S20]), its numerical estimation is θ ≈ 0.3.
(g-h) Average Lévy Lorentz gas (ALL).This model is emblematic for RWs with spatially-dependent steps; Its different properties are described in [S21, S22].We consider a RW on a 1d lattice with position dependent reflection or transmission probabilities r(k) or t(k).In the subdiffusive model (subALL), the transmission coefficient is taken to be In the continuous setting, this is equivalent to a space dependent diffusion coefficient for 0 . The persistent exponent is given by θ = 1 − 1/(3 + a) [S22].In the superdiffusive model (supALL), the reflection coefficient is taken to be In the continuous setting, this is equivalent to a space dependent diffusion coefficient for 0 . The persistent exponent is given by θ (i) Scaled Brownian motion (sBm).The sBm model represents RWs with time-dependent steps [S23-S26].Starting from X 0 t a process with i.i.d.symmetric jumps η 0 i with finite variance, we define the scale process of parameter β by X t ≡ X 0 ⌊t β ⌋ , or equivalently η t = A second way to define the sBm on discrete times is to consider steps η t following a binomial distribution of parameters (N t , 1/2) where N t is Poisson distributed of average λ(t) = t β−1 (the one used in Supplementary Figure 2).In the continuous setting, it amounts to a time-dependent overdamped Langevin equation, where η t is a white noise of unit variance.The persistence exponent of sBm is θ = β/2 [S24].
The sBm model allows one to compute the distribution of τ x0 ≡ T x0+∆x − T x0 analytically in the continuous limit, We plot this exact formula in Fig. 2 of the main text (x 0 = n, ∆x = 1) and in Supplementary Figure 3 (x 0 = m, ∆x = n).In the case of small ∆x, S26) in full agreement with the central result of this work, Eq. ( 1) of the main text.
B. Systematic numerical check of the scale-invariance of the time increments Based on the RW models described above, we confirm the scale invariance of the following random variables : • The number of records.We show in Supplementary Figure 2 that the number of records at time t is indeed scale invariant at large times as its average and standard deviation grow as expected as t 1/dw for all non-Markovian process considered.
• The time-increments between records.Supplementary Figure 3

C. Systematic numerical check of the asymptotic independence of record ages
Based on the RW models described above, we make systematic tests of the independence hypothesis between record ages: • In Supplementary Figure 4, we check the effective independence hypothesis used in Eq. (3) of the main text, by comparing P(max(τ m , . . ., τ m+n−1 ) ≤ T ) with n+m−1 k=m (1 − S(k, T )) for various values of n in the regime n dw and T small in comparison to m dw .The functional forms being the same for both distributions even for different values of n (by rescaling with n and m), this confirms numerically the approximation for all the non-Markovian models considered.
• In Supplementary Figure 5, we display the probability to make at least δ successive records (event called a record run of length δ) when n record runs have been performed.In other words, we look at the joint distribution {τ n = 1, . . ., τ n+k = 1}, which shows an exponential decay in the correlations between successive record ages, and thus provides an additional numerical check of the independence of τ k .We note that the time decay rate of the exponential varies with n for aging RWs, as expected from the dependency on the number n of records of the early time regime for the record age τ n . Supplementary

S3. DATA ANALYSIS
We provide a comprehensive description of the datasets used in this study, along with the general methodology that yields the results presented in the main text.Furthermore, we include complementary datasets that provide additional confirmation of our findings, including cases involving aging time series.
A. Details on the datasets used in the main text Here, we present the datasets used in the main text: (a) Elbe river discharge (m 3 /s).We consider the daily mean debit of the Elbe river measured in Dresden [S27].It was observed that this quantity presents correlation which can be modeled by subdiffusive non-Markovian RWs [S28].In this time series, we obtain H ≈ 0.14 by application of the DMA method [S29, S30], see below.
(b) Volcanic soil temperature ( o C).It was shown in [S31, S32] that the soil temperature monitored in the volcanic caldera of the Campi Flegrei area in Naples follows an fBm of parameter H ≈ 0.42 once the data have been detrended by removing the linear trends between two temperature extrema (between two solstices).Here we also detrend the data by removing the same linear seasonal trends.We display the data measured at the Monte Olibano (OLB) site.
(c) Trajectories of microspheres in agarose gel (nm).The trajectories [S33] represent the 2d motion of 50-nm polystyrene microspheres in agarose hydrogel (we consider that x and y displacement are i.i.d., giving us 2 independent 1d trajectories).There are 20 trajectories of 2000 frames which were analyzed in [S33] who obtained 1/d w ≈ 0.43.
(d) Motion of amoeba intracellular vacuoles (pixels = 106nm).We consider vacuole intracellular trajectories inside the amoeba in a 2d plane (we consider that x and y displacement are i.i.d., giving us 2 independent 1d trajectories) of at least 2048 frames.It was estimated in [S33] that the walk dimension verifies 1/d w ≈ 0.67.
(e) Trajectories of telomeres (µm).We use 2d trajectories [S33] of telomeres in the nucleus of untreated U2OS cells obtained in Ref. [S34] (we consider that x and y displacement are i.i.d., giving us 2 independent 1d trajectories).Similarly to Ref. [S33], we only consider trajectories where the mean-square displacement grows as t 0.5±0.05, which corresponds to 1/d w ≈ 0.25.
(f ) DNA RW on the Homo sapiens β-myosin heavy chain (HUMBMYH7).It was observed in [S35, S36] that the process is a RW with long-range correlations of Hurst exponent H ≈ 0.67.We estimate and remove the bias v = 1 N N t=1 η t in the data by replacing η t by η t − v. (g) Cumulative London air temperature ( o C.day).In this case, the temperature fluctuations (as for (b), we remove the linear trends between two solstices) are fractional Gaussian noise (fGn), in agreement to what was observed in [S37], of Hurst exponent H ≈ 0.8 obtained via the DMA method [S29, S30], see below.We consider the cumulative temperature fluctuation, which is then fBm.This quantity is of a particular interest in the studies of derivative pricing (see [S37]).
(h) Cumulative Ethernet traffic (10 bytes.ms).The dataset represents the number of packets going through an Ethernet cable every 10 ms at the Bellcore Morristown Research and Engineering facility [S38].In particular, in [S39], it was shown that the process is a fGn of dimension H ≈ 0.8.As for (e), the cumulative number of requests up to a given time t gives a non-Markovian process.We detrend the cumulative data as for (g) by removing the estimated bias v at every step.We display the measurement performed in August 1989.

B. Characterization and parametrization of the data used in the main text
In this section we provide the method developed to determine the walk dimension of the time series presented in the main text as well as numerical checks of their stationarity.
In order to obtain the walk dimension d w in a time series, one applies the celebrated Detrending Moving Average (DMA) [S29, S30] method, which consists in evaluating the typical fluctuations in a window of size ℓ regardless of any bias or deterministic trend.More precisely, for a dataset (X t ) t=0,...,N , we consider the windows of size

C. Analysis of complementary datasets
We additionally conducted an analysis of the following complementary datasets to further demonstrate the broad applicability of our results: (a) Volcanic soil temperature ( o C).Another dataset (measured at the Monte Sant'Angelo site) of soil temperatures monitored in the volcanic caldera of the Campi Flegrei area in Naples.It was shown in [S31, S32] that it is an fBm of parameter H ≈ 0.4 once the data have been detrended by removing the linear trends between two temperature extrema (between two solstices).Here we also detrend the data by removing the same linear seasonal trends.
(b) Cumulative air temperature at Montélimar ( o C.day).As for the soil temperature in [S31], we remove the linear trends between two solstices.In this case, the temperature fluctuations are fractional Gaussian noise (fGn) and not fBm, in agreement to what was observed in [S37] for the data at London.We obtain a Hurst exponent H ≈ 0.8 via the DMA method.
(c) Rhône river discharge (m 3 /s).We consider the daily mean debit of the Rhône measured at the Sault Brenaz station [S40].It was observed that this quantity presents correlations which can be modeled by subdiffusive non-Markovian RW [S28].In this time series, we obtain H ≈ 0.21 by applications of the DMA method.
(d) DNA RW.We consider the human T-cell receptor α/δ sequence from the GenBank data base (HUMTCRADCV).
It was observed in [S35, S36] that the process is a RW with long-range correlations of Hurst exponent H ≈ 0.61.We estimate and remove the bias v = 1 N N t=1 η t in the data by replacing η t by η t − v. (e) Cumulative Ethernet traffic (10 bytes.ms).The dataset represents the number of packets going through an Ethernet cable every 10 ms at the Bellcore Morristown Research and Engineering facility [S38].In particular, in [S39], it was shown that the process is a fGn of dimension H ≈ 0.84.This is why we consider the cumulative number of requests up to a given time t.We detrend the cumulative data in the same manner as the DNA RW by removing the estimated bias v at every step.We display the data measured in October 1989.
All considered complementary datasets support our theoretical results.